The following R packages will be needed:
library(sf)
library(dplyr)
library(tmap)
library(leaflet)
library(leafgl)
library(crosstalk)
library(sfdep)
library(ggplot2)
The homicide and geographic boundary data are loaded in. Note that 83 homicide data points had no geography associated with them and thus had to be dropped for spatial analysis. An additional 6 data points had no date, and had to be dropped as well.
# Load in homicide data and remove entries with missing geometry (n=83) and missing date (n=6)
homicides = read_sf("homicides_1940-1965_geocoded_v3/homicides_1940-1965_geocoded_v3.shp") %>%
filter(!st_is_empty(.)) %>%
filter(date != "NaT")
city_bounds = read_sf("Chicago_City_Limits-shp/Chicago_City_Limits.shp") %>%
st_transform(st_crs(homicides))
hexagons = read_sf("hexagons/hexagons.shp") %>%
st_transform(st_crs(homicides))
CAs = read_sf("Boundaries - Community Areas (current)/geo_export_e916a3b4-0ea3-4ca2-944a-cb060aacf705.shp") %>%
st_transform(st_crs(homicides))
A small amount of “jitter” is applied to points in order to make overlapping points both visible:
homicides_jitter = st_jitter(homicides, factor = 0.0001)
An interactive map of homicides:
leaflet(homicides_jitter) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(data = city_bounds, opacity = 0.2) %>%
addGlPoints(homicides_jitter, radius = 2, fillColor = "black")
Preliminary findings: Homicides 1940-1964 were concentrated primarily on the South and West Sides, with a smaller cluster appearing along the lake shore to the north as well. The Northwest and Southwest sides had a far lower concentration of homicides, though this may be due to lower population density at the time. On a more local scale, homicides tend to cluster around major thoroughfares, like Madison St on the West Side.
Spatial count of homicides per hexagon:
# Assign each point the grid ID of its surrounding hexagon
homicides_hex_joined = st_join(homicides, hexagons, join = st_within)
# Create table of frequency by hexagon
hex_freq_table = as.data.frame(table(homicides_hex_joined$GRID_ID)) %>%
rename(GRID_ID = Var1, hom_ct = Freq)
# Merge frequency table back into hexes
hexagons_hom_ct = hexagons %>%
left_join(hex_freq_table, by = "GRID_ID") %>%
mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))
Viewing the hexagons:
tmap_mode("view")
tm_shape(hexagons_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,Inf), title="Homicides", id="hom_ct", popup.vars = c("Grid ID:" = "GRID_ID", "Homicides:" = "hom_ct"))
Note that a custom classification scheme is used because of the extreme skew of the data. The 0 to 1 bin shows hexes where no homicide incidents occurred at all.
Spatial count of homicides per community area:
# Assign each point the grid ID of its surrounding CA
homicides_CA_joined = st_join(homicides, CAs, join = st_within)
# Create table of frequency by CA
CA_freq_table = as.data.frame(table(homicides_CA_joined$area_num_1)) %>%
rename(area_num_1 = Var1, hom_ct = Freq)
# Merge frequency table back into CAs
CAs_hom_ct = CAs %>%
left_join(CA_freq_table, by = "area_num_1") %>%
mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))
Viewing the community areas:
tmap_mode("view")
tm_shape(CAs_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,500,Inf), title="Homicides", id="hom_ct", popup.vars = c("Area Name:" = "community", "Homicides:" = "hom_ct"))
The same classification scheme is used for consistency, though it should be noted that community areas are generally larger than hexagons. Thus, homicides exceed 500 in certain areas, and so an extra bin, “500 or more”, is added here.
Preliminary findings: As seen above, homicides tend to cluster in the South and West Sides. However, while placement maps are generally good at highlighting areas of high concentration, the two choropleth maps above elucidate areas of low concentration as well, such as the community area of North Park, which had 0 homicides during the entire time period. These results should be interpreted with caution, however, as they are not population standardized, and may simply be apparent due to a lack of population density.
First, homicide dates must be converted into date format. Then
SharedData from crosstalk is used to link data between the
slider and the map below:
homicides_jitter$start = as.Date(homicides_jitter$date)
shared_data = SharedData$new(homicides_jitter)
The map:
filter_slider("date", "Date", shared_data, ~start, width = "100%")
leaflet(shared_data) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addPolygons(data = city_bounds, opacity = 0.2) %>%
addCircleMarkers(radius = 1, weight = 0.2, color = "black")
Preliminary findings: The interactive map above illustrates the outward spread of homicides in Chicago over time, likely paralleled by the city’s outward growth. However, the change in concentration of homicides in the West Side is quite notable; while in 1940 homicides appear to be no more concentrated on the West Side than any other area of the city, by 1964 they are just as concentrated in that area, if not more, than on the South Side.
Spatial count of homicides per hexagon for each year:
# Create table of frequency by hexagon and year
hex_freq_table_yr = as.data.frame(table(homicides_hex_joined$GRID_ID, homicides_hex_joined$year)) %>%
rename(GRID_ID = Var1, year = Var2, hom_ct = Freq)
# Convert factors to integers/strings
hex_freq_table_yr$GRID_ID = as.character(paste(hex_freq_table_yr$GRID_ID))
hex_freq_table_yr$year = as.integer(paste(hex_freq_table_yr$year))
# Create an expanded grid of hexagons that includes each hex/year combo
hexagons_yr = expand.grid(GRID_ID = hexagons$GRID_ID, year = 1940:1964) %>%
left_join(hexagons, by = "GRID_ID") %>%
st_as_sf()
# Merge frequency table back into hex/year grid
hexagons_hom_ct_yr = hexagons_yr %>%
left_join(hex_freq_table_yr, join_by(GRID_ID, year)) %>%
mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))
Creating a list of hexagon maps:
hex_map_list = list()
for (n in 1940:1964) {
data = hexagons_hom_ct_yr %>% filter(year == n)
map = tm_shape(data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, Inf), title = paste("Homicides in", n))
hex_map_list[[n-1939]] = map
rm(n, data, map)
}
Finally, the tmap_animation function is used to create
the animation presented below:
tmap_animation(hex_map_list, "hex_map_list.gif", width = 400, height = 500, delay = 50)
A similar set of code is used to create an animation of homicides per Community Area for each year:
# Create table of frequency by CA and year
CA_freq_table_yr = as.data.frame(table(homicides_CA_joined$area_num_1, homicides_CA_joined$year)) %>%
rename(area_num_1 = Var1, year = Var2, hom_ct = Freq)
# Convert factors to integers/strings
CA_freq_table_yr$area_num_1 = as.character(paste(CA_freq_table_yr$area_num_1))
CA_freq_table_yr$year = as.integer(paste(CA_freq_table_yr$year))
# Create an expanded grid of CAs that includes each CA/year combo
CAs_yr = expand.grid(area_num_1 = CAs$area_num_1, year = 1940:1964) %>%
left_join(CAs, by = "area_num_1") %>%
st_as_sf()
# Merge frequency table back into CA/year grid
CAs_hom_ct_yr = CAs_yr %>%
left_join(CA_freq_table_yr, join_by(area_num_1, year)) %>%
mutate(hom_ct = ifelse(is.na(hom_ct), 0, hom_ct))
Creating a list of community area maps:
CA_map_list = list()
for (n in 1940:1964) {
data = CAs_hom_ct_yr %>% filter(year == n)
map = tm_shape(data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste("Homicides in", n))
CA_map_list[[n-1939]] = map
rm(n, data, map)
}
Finally, the tmap_animation function is used to create
the animation presented below:
tmap_animation(CA_map_list, "CA_map_list.gif", width = 400, height = 500, delay = 50)
Preliminary findings: Because homicides are somewhat sporadic events over the course of a 1-year timespan, there appears to be a lot of oscillation between frames. This indicates that a longer timespan should be used for each frame. The animations below present the results with a 5-year timespan instead of 1-year.
hex_map_list_5yr = list()
CA_map_list_5yr = list()
for (n in seq(1940, 1964, 5)) {
data = hexagons_hom_ct_yr %>%
filter(between(year, n, n+4))
aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
rename(GRID_ID = Group.1, hom_ct = x) %>%
left_join(hexagons, by = "GRID_ID") %>%
st_as_sf()
map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste0("Homicides in ", n, "-", n+4))
hex_map_list_5yr[[1+((n-1940)/5)]] = map
rm(n, data, aggregated_data, map)
}
for (n in seq(1940, 1964, 5)) {
data = CAs_hom_ct_yr %>%
filter(between(year, n, n+4))
aggregated_data = aggregate(data$hom_ct, list(data$area_num_1), FUN = sum) %>%
rename(area_num_1 = Group.1, hom_ct = x) %>%
left_join(CAs, by = "area_num_1") %>%
st_as_sf()
map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, 100, Inf), title = paste0("Homicides in ", n, "-", n+4))
CA_map_list_5yr[[1+((n-1940)/5)]] = map
rm(n, data, aggregated_data, map)
}
tmap_animation(hex_map_list_5yr, "hex_map_list_5yr.gif", width = 400, height = 500, delay = 100)
tmap_animation(CA_map_list_5yr, "CA_map_list_5yr.gif", width = 400, height = 500, delay = 100)
Preliminary findings: Now the outward spread of homicides is a bit clearer, especially in the CA animation above. However, some oscillation is still present.
The following code creates and visualizes the standard deviational ellipse for the homicide data:
sde_data = std_dev_ellipse(homicides)
sde = st_ellipse(geometry = sde_data, sx = sde_data$sx, sy = sde_data$sy, rotation = -sde_data$theta)
leaflet(homicides_jitter) %>%
addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
addGlPoints(homicides_jitter, radius = 2, fillColor = "black") %>%
addPolygons(data = sde, opacity = 0.2)
Preliminary findings: Rather intuitively, the standard deviational ellipse encompasses the South and West Sides, with the Near North Side included as well. Future analysis could demonstrate the change of the ellipse’s size and shape over time, which could demonstrate how the geographic distribution of homicides spread out in a more systematic way.
A frequency table of homicides per hexagon before and after the time period shown above, compiled previously, is loaded in.
# Load in homicides per hexagon data
pt2_hex_freq_table_yr = read.csv("homicides_1874-2020_hex/homicides_pre_1931_post_1964.csv") %>%
select(GRID_ID, year, hom_ct = homicides)
First, all homicides 1874-1930 and 1965-2020 are visualized using a simple choropleth map. Note that records between 1874-1930 include all homicides; records used in the above analysis of homicides between 1940-1964 only include murders, justifiable homicides, and cases in which the verdict is unclear; and records between 1965-2016 only include murders. Thus, the findings from each of these time periods are not directly comparable to each other.
# Separate data between 1874-1930 and between 1965-2020
pt2_before_only = pt2_hex_freq_table_yr %>%
filter(year < 1931) %>%
select(GRID_ID, hom_ct)
pt2_after_only = pt2_hex_freq_table_yr %>%
filter(year > 1964) %>%
select(GRID_ID, hom_ct)
# Take aggregate of tables by GRID_ID
pt2_before_aggregated = aggregate(.~GRID_ID, data = pt2_before_only,FUN = sum)
pt2_after_aggregated = aggregate(.~GRID_ID, data = pt2_after_only, FUN = sum)
# Merge frequency table into hexes
pt2_before_hexagons_hom_ct = hexagons %>%
left_join(pt2_before_aggregated, by = "GRID_ID")
pt2_after_hexagons_hom_ct = hexagons %>%
left_join(pt2_after_aggregated, by = "GRID_ID")
Viewing the hexagons:
tmap_mode("view")
tm_shape(pt2_before_hexagons_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,Inf), title="Homicides 1874-1930", id="hom_ct", popup.vars = c("Grid ID:" = "GRID_ID", "Homicides:" = "hom_ct"))
tm_shape(pt2_after_hexagons_hom_ct) + tm_polygons(col = "hom_ct", breaks = c(0,1,5,10,50,100,Inf), title="Homicides 1965-2020", id="hom_ct", popup.vars = c("Grid ID:" = "GRID_ID", "Homicides:" = "hom_ct"))
The same custom classification scheme described above is used here for consistency, though without the “500 or more” bin because it is not used.
Preliminary findings: These two maps further demonstrate how homicides have spread outward over time, well beyond the time period of 1940-1964.
Next, time series animations like the ones shown above are created for hexagons between 1874-1930 and 1965-2016. A 5-year interval is used like before to reduce noise. Because of this, the years 1874, 1930, 2015, and 2016 will not be included, as they do not fit into any 5-year intervals.
# Create an expanded grid of hexagons that includes each hex/year combo
pt2_hexagons_yr = expand.grid(GRID_ID = hexagons$GRID_ID, year = c(1874:1930, 1965:2016)) %>%
left_join(hexagons, by = "GRID_ID") %>%
st_as_sf()
# Merge into hexagons
pt2_hexagons_hom_ct_yr = pt2_hexagons_yr %>%
left_join(pt2_hex_freq_table_yr, join_by(GRID_ID, year))
Creating a list of hexagon maps:
pt2_map_list_before = list()
pt2_map_list_after = list()
for (n in seq(1875, 1929, 5)) {
data = pt2_hexagons_hom_ct_yr %>%
filter(between(year, n, n+4))
aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
rename(GRID_ID = Group.1, hom_ct = x) %>%
left_join(hexagons, by = "GRID_ID") %>%
st_as_sf()
map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste0("Homicides in ", n, "-", n+4))
pt2_map_list_before[[1+((n-1875)/5)]] = map
rm(n, data, aggregated_data, map)
}
for (n in seq(1965, 2014, 5)) {
data = pt2_hexagons_hom_ct_yr %>%
filter(between(year, n, n+4))
aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
rename(GRID_ID = Group.1, hom_ct = x) %>%
left_join(hexagons, by = "GRID_ID") %>%
st_as_sf()
map = tm_shape(aggregated_data) + tm_polygons(col = "hom_ct", breaks = c(0, 1, 5, 10, 50, Inf), title = paste0("Homicides in ", n, "-", n+4))
pt2_map_list_after[[1+((n-1965)/5)]] = map
rm(n, data, aggregated_data, map)
}
Finally, the tmap_animation function is used to create
the two animations presented below:
tmap_animation(pt2_map_list_before, "hex_map_list_before.gif", width = 400, height = 500, delay = 100)
tmap_animation(pt2_map_list_after, "hex_map_list_after.gif", width = 400, height = 500, delay = 100)
Preliminary findings: The two animations above further demonstrate the spreading-out of homicides over time that parallels the outward growth of the city. The first animation, including homicides from 1865 to 1930, begins with homicides almost entirely concentrated in the city center. By 1925-1929, however, homicides are just as concentrated in the Near West and Near South Sides as they are in the city center. The latter animation also shows a diffusion of homicides over time, though also shows a decrease in their overall concentration within the last two time periods. It also suggests that the shape of this distribution has fundamentally changed over the past few decades – from a hub-and-spoke shape centered on the Loop and reaching into the South, West, and North Sides in the 1960s to two distinct regions on the West and South Sides and no central concentration in the Loop in the 2010s.
An animation of homicides per hexagon per 5-year interval for all time periods is created below.
tmap_animation(c(pt2_map_list_before, hex_map_list_5yr, pt2_map_list_after), "hex_map_list_all.gif", width = 400, height = 500, delay = 100)
Note: As mentioned above, methodologies varied between the three time periods (all homicides were included in 1874-1930; only murders, justifiable homicides, and unclear verdicts in 1940-1964; and only murders in 1965-2016) and thus time periods should be compared with caution.
In this section, the map of homicide counts per hexagon is reproduced, but with Zillow “ZTRAX” building counts used as a denominator, to get a rate of homicides based on how settled the area was at the time. An additional time series map animation of building counts is provided as well for reference.
hex_bldgs = read.csv("hexagons/hex_building_counts.csv")
pt2_hexagons_hom_ct_yr_ALL = rbind(hexagons_hom_ct_yr, pt2_hexagons_hom_ct_yr)
pt2_map_list_bldgs = list()
pt2_map_list_bldg_weighted = list()
for (n in seq(1875, 2014, 5)) {
if (n < 1930 | n > 1939) {
# Homicide data
data = pt2_hexagons_hom_ct_yr_ALL %>%
filter(between(year, n, n+4))
aggregated_data = aggregate(data$hom_ct, list(data$GRID_ID), FUN = sum) %>%
rename(GRID_ID = Group.1, hom_ct = x) %>%
left_join(hexagons, by = "GRID_ID") %>%
st_as_sf()
# Building data
bldgs = hex_bldgs %>%
filter(between(year, n, n+4))
aggregated_bldgs = aggregate(bldgs$weighted_building_count, list(bldgs$GRID_ID), FUN = mean) %>%
rename(GRID_ID = Group.1, bldg_ct = x)
# Combined data
aggregated_data = aggregated_data %>%
left_join(aggregated_bldgs, by = "GRID_ID") %>%
mutate(hom_rt = hom_ct / bldg_ct) %>%
mutate(hom_rt = ifelse(hom_rt == Inf, NA, hom_rt))
# Maps
bldgs_map = tm_shape(aggregated_data) + tm_polygons(col = "bldg_ct", breaks = c(0,1,5,10,50,100,500,1000,Inf), title = paste0("Buildings in ", n, "-", n+4))
map = tm_shape(aggregated_data) + tm_polygons(col = "hom_rt", breaks = c(0, 0.001, 0.01, 0.1, 1, 10), title = paste0("Homicides per building in ", n, "-", n+4))
pt2_map_list_bldgs[[1+length(pt2_map_list_bldgs)]] = bldgs_map
pt2_map_list_bldg_weighted[[1+length(pt2_map_list_bldg_weighted)]] = map
rm(n, data, aggregated_data, bldgs, aggregated_bldgs, map, bldgs_map)
}
}
An animation of building-weighted homicide rates per hexagon per 5-year interval is presented below:
tmap_animation(pt2_map_list_bldg_weighted, "hex_map_list_bldg_weighted.gif", width = 400, height = 500, delay = 100)
Note: Areas labeled “missing” had no buildings, so the calculation of the rate produced a division by 0 error.
tmap_animation(pt2_map_list_bldgs, "hex_map_list_bldgs.gif", width = 400, height = 500, delay = 100)
In this section, time series line graphs are made to show the rate of homicides per building for each year in hexagons categorized by annexation status and land sale year.
Previous homicide and building data are used to make a dataset for Part 3 that characterizes each hexagon by its homicide count, building count, and homicide per building rate for each year.
pt4_data = left_join(pt2_hexagons_hom_ct_yr_ALL, hex_bldgs, by = c("GRID_ID", "year")) %>%
select(GRID_ID, year, hom_ct, bldg_ct = weighted_building_count) %>%
mutate(hom_rt = hom_ct / bldg_ct) %>%
mutate(hom_rt = ifelse(hom_rt == Inf, NA, hom_rt))
First, hexagons are categorized by the following classification scheme: * Hexagons within the original Chicago city limits * Hexagons within areas incorporated into the city by the state * Hexagons within settlements that were annexed into the city * Hexagons within settlements (Jefferson and Lake Townships) that were sparsely populated when they were annexed into the city.
This is done by calculating the centroids of the hexagons, and spatially assigning those centroids to the surrounding area. Additionally, any hexagon centroids that are outside of the city limits are excluded.
annexed_areas = read_sf("annexed/annexed_areas/Annexed Settlements.shp") %>%
st_transform(st_crs(hexagons))
hexagons_annexed = hexagons %>%
st_centroid() %>%
st_filter(city_bounds, .predicate = st_within) %>%
st_join(annexed_areas, join = st_within) %>%
select(-year) %>%
mutate(annex_category = case_when(
settlement == "CHICAGO" ~ "Original",
is.na(settlement) ~ "Incorporated",
settlement == "LAKE" | settlement == "JEFFERSON" ~ "Annexed - Sparse",
.default = "Annexed - Populated"
)) %>%
st_drop_geometry() %>%
right_join(hexagons, ., by="GRID_ID") %>%
select(GRID_ID, settlement, annex_category)
The map below summarizes the categorization of hexagons by annexation status:
tmap_mode("view")
tm_shape(hexagons_annexed) + tm_polygons(col = "annex_category", title="Annexation Category", id="annex_category", popup.vars = c("Grid ID:" = "GRID_ID", "Settlement:" = "settlement", "Annexation Category:" = "annex_category"))
Next, the annexation table is joined to the Part 3 dataset to categorize each observation by annexation status. Then, observations are aggregated by annexation status to produce the time series graph shown below.
pt4_annexed = right_join(pt4_data, st_drop_geometry(hexagons_annexed), by = "GRID_ID") %>%
st_drop_geometry()
pt4_annexed_aggregated = pt4_annexed %>%
group_by(year, annex_category) %>%
summarise(hom_ct = sum(hom_ct), bldg_ct = sum(bldg_ct)) %>%
mutate(hom_rt = hom_ct / bldg_ct)
# Line Graph
ggplot(data = pt4_annexed_aggregated, aes(x = year, y = hom_rt, group = annex_category)) +
geom_line(aes(color = annex_category), size = 1) +
ylab("Homicides per Building") +
xlab("Year") +
coord_cartesian(ylim = c(0, 0.01))
Note: The y-axis is limited at 0.01 to reduce the impact of noise on the interpretability of the rest of the graph.
The animation below shows homicide rates by year, separated by annexation category.
mean_na_rm = function(x) {
mean(x, na.rm = TRUE)
}
pt4_map_list_annexed = list()
for (n in seq(1875, 2014, 5)) {
if (n < 1930 | n > 1939) {
data = pt4_annexed %>%
filter(between(year, n, n+4))
aggregated_data = aggregate(data$hom_rt, list(data$GRID_ID), FUN = mean_na_rm) %>%
rename(GRID_ID = Group.1, hom_rt = x) %>%
left_join(hexagons_annexed, by = "GRID_ID") %>%
st_as_sf()
map = tm_shape(aggregated_data) + tm_facets(by = "annex_category", free.coords = F, scale.factor = 2) + tm_polygons(col = "hom_rt", breaks = c(0, 0.001, 0.01, 0.1, 1, 10), title = paste0(n, "-", n+4))
pt4_map_list_annexed[[1+length(pt4_map_list_annexed)]] = map
rm(n, data, aggregated_data, map)
}
}
An animation of building-weighted homicide rates per hexagon per 5-year interval is presented below:
tmap_animation(pt4_map_list_annexed, "hex_map_list_annexed.gif", width = 800, height = 1000, delay = 100)
Work in progress
Finally, hexagons are categorized by redlining grades, A, B, C, D, and other (includes industrial, commercial, or sparsely built-up areas).
This is done by calculating the centroids of the hexagons, and spatially assigning those centroids to the surrounding area’s redlining grade. Additionally, any hexagon centroids that are outside of the city limits are excluded.
redlined_areas = read_sf("redlining/redlining.geojson") %>%
st_transform(st_crs(hexagons)) %>%
st_make_valid() %>%
select(grade)
hexagons_redlined = hexagons %>%
st_centroid() %>%
st_filter(city_bounds, .predicate = st_within) %>%
st_join(redlined_areas, join = st_within) %>%
st_drop_geometry() %>%
right_join(hexagons, ., by="GRID_ID") %>%
select(GRID_ID, grade)
The map below summarizes the categorization of hexagons by annexation status:
tmap_mode("view")
tm_shape(hexagons_redlined) + tm_polygons(col = "grade", palette = c("#00aa00", "#0000ff", "#ffff00", "#ff0000"), title="Redlining Grade", id="grade", popup.vars = c("Grid ID:" = "GRID_ID", "Redlining Grade:" = "grade"))
Next, the redlining table is joined to the Part 3 dataset to categorize each observation by redlining grade. Then, observations are aggregated by redlining grade to produce the time series graph shown below.
pt4_redlined = right_join(pt4_data, st_drop_geometry(hexagons_redlined), by = "GRID_ID") %>%
st_drop_geometry()
pt4_redlined_aggregated = pt4_redlined %>%
group_by(year, grade) %>%
summarise(hom_ct = sum(hom_ct), bldg_ct = sum(bldg_ct)) %>%
mutate(hom_rt = hom_ct / bldg_ct)
# Line Graph
ggplot(data = pt4_redlined_aggregated, aes(x = year, y = hom_rt, group = grade)) +
geom_line(aes(color = grade), size = 1, palette = c("#00aa00", "#0000ff", "#ffff00", "#ff0000")) +
ylab("Homicides per Building") +
xlab("Year") +
scale_color_manual(values = c("#00aa00", "#0000ff", "#ffff00", "#ff0000")) +
coord_cartesian(ylim = c(0, 0.01))
Note: The y-axis is limited at 0.01 to reduce the impact of noise on the interpretability of the rest of the graph.
The animation below shows homicide rates by year, separated by annexation category.
pt4_map_list_redlined = list()
for (n in seq(1875, 2014, 5)) {
if (n < 1930 | n > 1939) {
data = pt4_redlined %>%
filter(between(year, n, n+4))
aggregated_data = aggregate(data$hom_rt, list(data$GRID_ID), FUN = mean_na_rm) %>%
rename(GRID_ID = Group.1, hom_rt = x) %>%
left_join(hexagons_redlined, by = "GRID_ID") %>%
st_as_sf()
map = tm_shape(aggregated_data) + tm_facets(by = "grade", free.coords = F, scale.factor = 2) + tm_polygons(col = "hom_rt", breaks = c(0, 0.001, 0.01, 0.1, 1, 10), title = paste0(n, "-", n+4))
pt4_map_list_redlined[[1+length(pt4_map_list_redlined)]] = map
rm(n, data, aggregated_data, map)
}
}
An animation of building-weighted homicide rates per hexagon per 5-year interval is presented below:
tmap_animation(pt4_map_list_redlined, "hex_map_list_redlined.gif", width = 800, height = 1500, delay = 100)
Work in progress